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Abstract 

In this paper, we study the numerical approximation of a system of par¬ 
tial differential equations describing the corrosion of an iron based alloy 
in a nuclear waste repository. In particular, we are interested in the con¬ 
vergence of a numerical scheme consisting in an implicit Euler scheme in 
time and a Scharfetter-Gummel finite volume scheme in space. 
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1 Introduction 

1.1 General framework of the study 

At the request of the French nuclear waste management agency ANDRA, investi¬ 
gations are conducted to evaluate the long-term safety assessment of the geological 
repository of high-level radioactive waste. The concept of the storage under study 
in France is the following: the waste is confined in a glass matrix, placed into cylin¬ 
drical steel canisters and stored in a claystone layer at a depth of several hundred 
of meters. The long-term safety assessment of the geological repository has to take 
into account the degradation of the carbon steel used for waste overpacks, which is 
mainly caused by generalized corrosion processes. 

In this framework, the Diffusion Poisson Coupled Model (DPCM) has been pro¬ 
posed by C. Bataillon et al. in [3] in order to describe corrosion processes at the 
surface of the steel canisters. It assumes that the metal is covered by a dense oxide 
layer which is in contact with the claystone. The model describes the evolution of 
the dense oxide layer. In most industrial cases, the shape of metal pieces is not flat 
(container or pipe). But the oxide layer is very thin compared to the size of the 
exposed surface. In practice, the available data are averaged over the whole exposed 
surface and the microscopic or even macroscopic heterogeneities of materials are not 
taken into account. For these reasons, a ID modeling has been proposed in order to 
describe the real system. 

The oxide layer behaves as a semiconductor: charge carriers, like electrons, 
cations Fe 3+ and oxygen vacancies, are convected by the electric held and the elec¬ 
tric potential is coupled to the charge densities. The DPCM model is then made 
of drift-diffusion equations on the charge densities coupled with a Poisson equation 
on the electric potential. The boundary conditions induced by the electrochemical 
reactions at the interfaces are Robin boundary conditions. Moreover, the system 
includes moving boundary equations. 

Numerical methods for the approximation of the DPCM model have been de¬ 
signed and studied by Bataillon et al in pQ. Numerical experiments with real-life 
data shows the ability of the model to reproduce expected physical behaviors. How¬ 
ever, the proof of convergence of the scheme is challenging. In this paper, we will 
focus on a simplified model with only two species: electrons and cations Fe 3+ . As 
the displacement of the interfaces in the full DPCM model is due to the current 
of oxygen vacancies, which are not taken into account in this study, the simplified 
model will be posed on a fixed domain. These simplifications will permit us to show 
how to deal with the boundary conditions for the numerical analysis of the scheme 
introduced in [Tj. 

1.2 Presentation of the model 

In this paper we focus on the simplified corrosion model already introduced in mm- 
The unknowns of the model are the densities of electrons N and cations Fe 3+ P and 
the electric potential T. The current densities are respectively denoted by Jjy and 
Jp; they contain both a drift and a diffusion part. Therefore, the model consists in 
two drift-diffusion equations on the densities, coupled with a Poisson equation on 
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the electric potential. Let T > 0, the model is written as 

d t P + d x J P = 0, J P =-d x P - 3Pd x ^>, in (0,1) x (0,T), (la) 
ed t N + d x J N = 0, J N = -d x N + Nd x ^>, in (0,1) x (0, T), (lb) 

-\ 2 d% x 'f! = 3P - N + p h i, in (0,1) X (0,T), (lc) 

where A is the rescaled Debye length and phi the net charge density of the ionic 
species in the host lattice which is constant. The parameter e stands for the ratio 
of the mobility coefficients of electrons and cations, then e -C 1. 

As the equations ()la[) and (llbl) for charge carriers densities have the same form, 
we will use the following synthetical form: 

e u d t u + d x J u = 0, J u = —d x u - z u ud x T, in (0,1) x (0, T). (2) 

For u = P,N, the charge numbers of the carriers are respectively z u = 3, —1 and we 
respectively have e u = l,e. 

Let us now focus on the boundary conditions. Charge carriers are created and 
consumed at both interfaces x = 0 and x = 1. The kinetics of the electrochemical 
reactions at the interfaces are described by Butler-Volrner laws. It leads to Robin 
boundary conditions on N and P. As in [T], we assume that the boundary conditions 
for P and N have exactly the same form. Therefore, for u = P,N, they are written 
as 


—J u = r®(u, T) on {x = 0} x (0, T), (3a) 

J u = rl(u, T, V) on {x = 1} x (0,T), (3b) 

where V is a given applied potential (we just consider here the potentiostatic case) 
and and r\ are linear and monotonically increasing functions with respect to their 
first argument. More precisely, due to the electrochemical reactions at the interfaces, 
we have, for u = P,N: 


rl(s,x) = ffi(x)s- 7 °(x), (4a) 

ri(a,x,V) = $(V-x)s-'yi(V-x), (4b) 


where the functions (j3 l u )i= o,i, (7„)j=o,i are given functions. These functions depend 
on many parameters: the interface kinetic coefficients (m^, ^)j=o i, the positive 
transfer coefficients (a l u , 14 )*=o,i, the maximum occupancy for octahedral cations in 
the lattice p max and the electron density in the state of metal N max . For u = P,N, 
the functions (/3„)j=o,i, (7«)i=o,i are written 


= m l u e 


—Zublx 


+ Ke z 


* = 0 , 1 , 


7 °( x ) = rn 0 u u max e - Zub ° x , 7 i(x) = kiu max e z -< x . 


(5a) 

(5b) 


Throughout the paper, we will assume that the interface kinetic coefficients and 
the transfer coefficients are given constants which satisfy 


m. 


°ik° u ,m l u ,kl> 0, for u = P,N, 


( 6 ) 
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a ui b h a l, b l € [0,1], for u = P,N. (7) 

We also assume that phi does not depend on x and that 

3 pmax _ N max + ^ = ( 8 ) 

Indeed, in the applications (see 0), the scaling of the model leads to phi = —5, 
pmax _ 2 anc [ jymax _ ^ go ^hat ^he relation © is satisfied. 

The boundary conditions for the Poisson equation take into account that the 
metal and the solution can be charged because they are respectively electronic and 
ionic conductors. Such an accumulation of charges induces a field given by the Gauss 
law. These accumulations of charges depend on the voltage drop at the interface 
given by the usual Helmholtz law which links the charge to the voltage drop through 
a capacitance. The parameters ATq~ c and AT^' C are the voltage drop corresponding 
to no accumulation of charges respectively in the metal and in the solution. Finally, 


the boundary conditions for the electric potential are written: 

T — aod x 'k = A'Pq ZC , on {x = 0} x (0, T), (9a) 

T + aid x V = V- ATf c , on {x = 1} x (0,T), (9b) 

where ao and aq are positive dimensionless parameters arising in the scaling. 

The system is supplemented with initial conditions, given in L°°(0,1): 

u(x, 0) = u°(x), for u = P,N. (10) 

Moreover, we assume that these initial conditions satisfy 

0 ^ u° < u max , a.e. on (0,1), for u = P,N. (11) 


In the remainder of the paper we will denote by (V) the corrosion model defined 
by ©, m, ©, m and m■ Let us first note that the system of equations © 
is the so-called linear drift-diffusion system. This model is currently used in the 
framework of semiconductors device modeling (see for instance [371 ESI EH [33]) or 
plasma physics (see [15]). In this context the drift-diffusion model has been widely 
studied, from the analytical as from the numerical point of view. Let us refer 
to the pioneering work by Gajewski [20] about existence and uniqueness results. 
Further developments have been done in mmmm- The long-time behavior 
of solutions to the drift-diffusion system via an entropy method has been studied in 
[22 , 27] and the stability at the quasi-neutral limit or at the zero-electron-mass limit 
in [24. 1291 30, 31]. Different methods have also been proposed for the approximation 
of the drift-diffusion system, and studied, see for instance mm for finite volume 
schemes, [35] for a mixed finite volume scheme and [SI E El EES] for finite element 
and mixed exponential fitting schemes. 

In the modeling of semiconductor devices, the boundary conditions are generally 
mixed Dirichlet/Neumann boundary conditions (corresponding to the ohmic con¬ 
tacts and the insulated boundary segments of the device). Then, the originality of 
the corrosion model described in this paper lies in the boundary conditions ©, © 
which are of Robin type, and induce an additional coupling between the equations. 
Let us define the notion of weak solution to the corrosion model ( V ). 
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Definition 1.1 We say that (P, N, T) € L 2 (0, T; i7 1 (0,1)) n L°°([0,T] x [0,1]) is a 
weak solution of (V), if for all ip in L 2 (0, T; i/ 1 (0,1)), 


— £, 


/*T rl i* 1 rT rl 

, / / udtpdxdt — e u / uo(x)</?(0, x)dx — / / (—d x u — z u ud x ^) d x <pdxdt 

Jo Jo Jo Jo Jo 

+ / [(/^(^ - !))«(*, !) - 7 «(^ - 1 ))) ¥>(*, 1 ) 

Jo 

+ - 7°O0,0))) 0)] dt = 0, for u = P,N, (12) 


and 

A 2 [ f d x 'S?d x (pdxdt — [ — (V — — A^ zc ) ip(t,l)dt 

Jo Jo Jo a i 

+ f — (\F(t, 0) - A\E^~ C ) ip(t, 0)dt = [ f (3P - N + p h i)pdxdt. (13) 
Jo a o Jo Jo 

In ED, Chainais-Hillairet and Lacroix-Violet have proved, under some assump¬ 
tions on the chemical and physical parameters, the existence of a weak solution to 
(' V). This result is obtained by passing to the limit in an approximate solution given 
by a semi-discretization in time. Convergence of the sequence of approximate solu¬ 
tions is ensured by some estimates which yield compactness. In the current article, 
we apply the same ideas as in ED to a full discretization of (V) presented below. 

1.3 Presentation of the numerical scheme 

Some numerical schemes for the approximation of ( V) have been proposed by Batail- 
lon et al in [T], Their stability analysis is fulfilled but no convergence results are 
given. Here we are interested in the convergence analysis of the fully implicit scheme 
introduced in [lj. It is a backward Euler scheme in time and a finite volume scheme 
in space, with Scharfetter-Gummel approximation of the convection-diffusion fluxes. 

Let us consider a mesh T for the domain [0,1]. It consists in a family of mesh 
cells denoted by for i € [l;/], with 

0 = x l/2 < x 3/2 < ■ ■ ■ < Xj_ 1/2 < x I+1/2 = 1. 

Then, we define x % = l+1 / 2 — 1 1 / 2 ; f or { g [1;/] and xq = x^/i = 0, x /+1 = 
x 1 + 1/2 = 1- Moreover, we set 

h i = x i+ i-x i _i, [l;/], 

h i+ i = x i+1 - x i} Vie [0; /]. 

The mesh size is defined by h = max {hi, i € [1: /]}. 

Let us denote by At the time step. We will assume that there exists I\ € N such 
that KAt = T (either, we would define K as the integer part of T/At). We consider 
the sequence {t k )o^k^K such that t k = kAt. 
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The scheme under study in this paper is written as follows. For i E [[1;/], 
A'e 10;/<- lj. 


—A z ( <W k +l - dT ft +( ) = hi (3ftf +i - Af +1 + p hl 


i+i 


yk+1 

' i 


jfc + l ATk+1 


? .^+l _ q.k 

e u hi r-i- + ft fc+1 ! - P k+l ! = 0, for u = P , AT, 

At u ^+\ 


with the numerical fluxes defined for i E [0; /]] by 


dT^+j 1 = 


' i+1 


h, 


*+5 


Trfc+1 


B (*A + .d**«) u, t+1 - B 


'■h-4 


(14a) 

(14b) 

(15a) 

for u = P, IV, 
(15b) 


where 1? is the Bernoulli function : 

B(x) = ———, Vx / 0 and 11(0) = 1. 
e® — 1 

We supplement the scheme with the discretization of the boundary conditions: for 


k € [0; ft' — 1], 

d> k+1 - cxodTi" 1 " 1 = AT^ C , (16a) 

2 

T^+i 1 + aid^+i = V - ATf c , (16b) 

J '2 

-^V = ^ (^o +1 ) *4 +1 - 7^ (^o +1 ) for u = P,N, (16c) 

= A (V ~ u k X\ — 7u (y — ^+!) for u = P,N, (16d) 
and of the initial conditions: for i E [[1; ft]], 

u i = [ + ^u°(x)dx, for u = P,N. (17) 

hi Jx. 1 
*“2 

The scheme (I14I) - (I17|) will be denoted in what follows by ( S). 


Remark 1 The choice of the Bernoulli function for B corresponds to a Scharfetter- 
Gummel approximation of the convection-diffusion fluxes. These numerical fluxes 
have been introduced by Il’in in [25] and Scharfetter and Gummel in [36] for the 
numerical approximation of the drift-diffusion system arising in semiconductor mod¬ 
elling. Lazarov, Mishev and Vassilevsky in [32] have established that they are 
second-order accurate in space. Dissipativity of the Scharfetter-Gummel scheme 
with a backward Euler time discretization for the classical drift-diffusion system was 
proved in [22] and Chatard in [l4j . One crucial property of the Scharfetter-Gummel 
fluxes is that they generally preserve steady-states. 
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Remark 2 The scheme (5) is written on a nonuniform discretization of the domain 
[0,1]. Indeed, the numerical experiments done in pLOj for the steady-state of ( S ) 
show some boundary layers for the density profiles. Therefore, it seems relevant to 
use a mesh which is refined near the boundaries. We will use a Tchebychev mesh, 
already introduced in [ID] and jlj. 

For the discretization in time, it is easier to deal with a fixed time step when 
studying the convergence. However, it would be also possible to introduce a variable 
time step. Adaptive time step strategy could also be taken into account, see [lj. 

1.4 Main results 

The aim of this paper is to prove the convergence of a sequence of approximate 
solutions obtained with the numerical scheme ( S ) to a solution of ( V ). 

To this end, we first need to establish the existence of a solution to the scheme. 
Indeed, at each time step k € [0;AT — 1], the vector of discrete unknowns 
( P fc+ \ N fc+ \4> fc+ i), with P k+1 = (P k+1 ) o^j+i, N k+1 = (N k+1 ) o^j+i , 
\F fc+1 _ i s defined as a solution to the nonlinear system of equa¬ 

tions (EMEU. In m, the existence of a solution has been proved only in the case 
where e > 0. As stated in ProDosition ll.il the result holds even if e = 0. 

Proposition 1.1 Let e > 0, ao > 0, a\ > 0 and the hypotheses ©, (0), ([8jh (fill) 
hold. Let us assume that 

- xV i 1 + lo § (ao apk° P )) ^ A4^ c < -jj- (! + log ( a 0 a° N k ° N )) , (18) 

oCLp d'N 

~ -r (1 + log (ai^mjv)) sj ATf c < i 1 + lo § (“i b P m p)) • ( 19 ) 

Then there exists a solution (P fc+1 , N fc+1 . 4 ,fc+1 )o<fc<A'-i to the fully implicit scheme 
(5). Moreover, it satisfies the following stability property: 

0 < if ^ P max and 0 < N k SC N max , Vi G [1; I + 1], VA; G |0; K\. (20) 

Based on the vector of discrete unknowns, we can dehne some approximate 
solutions, which are piecewise constant function in space and time, as it is usual 
for finite volume approximations. For a given mesh T of size h and a given At, we 
dehne, for w = N,P or T, 

i 

w h = y " 1 ! r | + W)H{x=o} + w k +1 t{ x= i}, for k € [0; A], 

i= 1 
K -1 

WhAt = w h +lt [t*p +!)• 
k =0 

For a sequence of meshes and time steps (T m ,A f m ) m such that h rn — >• 0 and 
A t m —» 0 as m —>• +oo, we can dehne a sequence of approximate solutions 
(Am, Nmi ^ ra)m with w m = Wh m Atm I° r w = P,N or T. The main result of the 
paper is the convergence of such a sequence of approximate solutions to a weak 
solution of the corrosion model ( V ). It is given in Theorem ll.il 


( 21 ) 

( 22 ) 
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Theorem 1.1 Let e > 0, «o > 0, a\ > 0. Assuming ©, (0), (|8|). (fill) . (fl8l) and 
m, there exist P, N and T £ L 2 (0,T; H 1 ( 0,1)) such that, up to a subsequence, as 
m —>• Too, 


P m -1 P strongly in L 2 (0, T; L 2 (0, 1)), 
Nm^-N strongly in L 2 (0,T;L 2 (0,1)), 
^ T strongly in L 2 (0,T;L 2 (0,1)), 


where (P, N, T) is a weak solution of (V) in the sense of Definition ll.il 

As it is well known in the finite volume framework (see for instance US]), the proof 
of Theorem o will be based on estimates satisfied by the approximate solutions 
and on compactness results. Due to the particular boundary conditions, we also 
need some additional results on the convergence of traces. They will be obtained 
following the ideas of [6]. 

The paper is organized as follows. Section 2 is devoted to the proof of Proposi¬ 
tion [lj] and also to the proof of discrete L 2 (0, T, if ^-estimates on the approximate 
densities and on the approximate potential. Then, in Section 3, we establish the 
compactness of the sequences of approximate solutions. We also obtain a conver¬ 
gence result for the traces on the boundaries. In Section 4 we conclude the proof of 
Theorem 11.11 by passing to the limit in the numerical scheme. Finally, Section 5 is 
devoted to the presentation of some numerical experiments. 


2 Existence result and discrete L 2 (0, T, i/ 1 )-estimates 


2.1 Proof of Proposition 11.11 

The proof of Proposition 11.11 for e > 0 has already been done in |1J. Then, we only 
prove the result for e = 0. To this end, we follow the ideas of [1] and [5]. 

Letting p > 0, we introduce a mapping Tjf : R /+2 x R /+2 —>• R /+2 x R /+2 such 
that 7 ~i^(P, N) = ( P , N ). This mapping is based on a linearization of the scheme ; it 
is defined in two successive step. First, we compute T as the solution to the linear 
system 

-A 2 (dtf i+ i - d,_x) = hi (3 P -Ni + p h i ), Vi € [1; /], 

with a definition of the numerical fluxes dT ,. + 1 analog to (I15al) and boundary condi¬ 
tions similar to (I16al) - (ll6bl) . The matrix of the linear system is obviously invertible 
and T is uniquely defined. 

Then, we define P and N as the solution to the linear systems 


hj 

At 


1 + 


Pi-^Pi-pn+p, 


A 2 


P, i+- 


- w 


P,i 


- 1 = 0 , 

2 


Vi€ [1;/], 


hi p_ 
At A 2 



+ i+1 


Vie [1;/], 
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with a definition of the numerical fluxes T u i+ i analog to (|15bl) and boundary con¬ 
ditions similar to (I16cl) - ()16dl) . As shown for instance in [I], the matrices defined at 
this step are M-matrices. Therefore, they are invertible and, as in [Tj, we can deduce 
that Tjf preserves the set 

K = { (P, N) G R /+2 x R /+2 ; 0 ^ Pi ^ P max , 0 ^ N t ^ N max , VO ^ i ^ I + l} , 

as long as (fT8l) - (fT9l) are satisfied and At verifies: 

At ^ umin (-, - ] . (23) 

V QP max Mmax I v ' 

Finally, is a continuous mapping from R /+2 x R^ +2 to itself which preserves 
the set 1C. Thanks to Brouwer’s Theorem, we conclude that 7)f has a fixed point in 
1C. This fixed point with the corresponding T defines a solution to (5) with e = 0 
at time step k + 1. Since ji is an arbitrary constant, we can choose it such that (1231) 
is verified and a solution to the scheme (5) satisfies (12011 without any condition on 
At. 


2.2 Notations and preliminary results 

In order to prove Theorem 11.11 we need to define some functional sets and norms 
and to establish some properties. This is the goal of this section. 

First of all, let us define two sets of functions. 

Definition 2.1 Let T be a mesh of size h and At a discretization time step. We 
first define Tij- the set of piecewise constant functions in space as 

7ir = {wh ■ [0,1] R | 3(tt!;)o<j</+i € R /+2 such that 

i 

u>h(x) = X) u,il (s i _ 1/2 ,*i +1/2 )(®) +™ O l{x=o}0*0 + w I+ it {x=l} (x) 

2=1 

Then, we define the set of piecewise constant functions in space and time 

y.T,At = \w h ,At '■ [O' !] X [0,T] H- R I 3(wJ; +1 )o<fc<A-_i G (' Ur) K such t! 

K—l 

Wh,At(x,t) = w h + 1 ( x ) 1 [t k ,t k + 1 )(. t ) 

k =0 



Then, denoting by || • ||o the usual L 2 (0, l)-norm, we remark that 


i \ V2 

..2 


IKIlo = Y hiw'i Vwh G Ut- 


\i=l 
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Moreover, we define some norms on 'Hj- and Ht, At, which are discrete counterparts 
of H^O,1), H~ l { 0,1) and L 2 (0, T, l))-norms: 


lkft,||i,T = + w 2 i+i] e Hr, 

\ 2=0 *+| ) 

II Wh II - 1 , 2 , r = max jy whVhdx, v h G Hr and ||r7i||i,T ^ l| € Hr, 

||Wfc,Ai||o ; l,T = f ^ ^ i ll tl; ft + 1 |ll,r') V w h,At € ^T, At- 


\fc =0 / 

As shown in [2j, for all wr € Hr, we have: 

(wi ) 2 < 211^11^7- Vi € [1; FJ, (26) 

and, as a direct consequence, the following discrete Poincare inequalities: 

IKIIo^V2|Klli,r Vur e H r - (27) 


Finally, we end this section with some properties satisfied by the functional 
spaces. These properties will be crucial in order to apply some compactness results. 
More precisely, the following lemmas state that the hypotheses (HI) and (H2) of 
Lemma 3.1 in [23] hold for sequences (Hr™)™ and (wh m )m such that for all m, 
w h m e Hr m . 


Lemma 2.1 Let {'H% n )m. be a sequence of finite-dimensional subspaces of L 2 (0,1) 
defined by ([Ml) . Let (ur m ) m be a sequence such that Wh m € U-Tm for all m and 
satisfying : 

3C > 0 such that Vm, \\w hrri \\i,T m ^ C. 

Then, up to a subsequence, (■ Wh m ) m converges to in L 2 (0,1) when m tends to 
Too. 


The proof of this lemma is a consequence of a Kolmogorov’s compactness Theorem 
(see for instance Theorem 10.3 in [18]). 


Lemma 2.2 Let (Hr^m be a sequence of finite-dimensional subspaces of L 2 (0,1) 
defined by (|24|) . Let (wh m )m be a sequence such that Wh m € W-Tm for all m. If 
( w hm)m converges to w in L 2 (0,1) and \\wh m \\-i,2,Tm converges to 0, then w = 0. 


Proof We will obtain that w = 0, as a consequence of: 

\hp € C“([0,1]), f w(f dx = 0. 

Jo 

Thus, letting ip G C“([0,1]), we set ifi = <p(xi) for all i G flO; I T lj and we define 
the associate ip hm G H Tm - Since w hm and <p hrn G H Tm : 



W hm^Phm 


dx 


^ \\ W h 7n \\-l,2 ) Tm\\^h m \\l ) Tm.- 


(28) 
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But, thanks to the regularity of ip, 


Wh m vdx 


€ 


[ W hm < Ph m dx 

Jo 


+ 


w h m (<P~ Vh m )dx 


^ 11 ^Xh m || 1,2, 'Tin 11 hm ||l, 7m + 11 w h m 11 oC'iphm ■ 

Since ||o is bounded independently of h m , for all p E C£°([0,1]), we have 


/ 


unpdx = hm / Wh m (pdx = 0, 


m—>•+oo 


and then w = 0. 


□ 


2.3 Discrete L 2 (0, T, i/ 1 (0,1)) estimates on P, N and 'k 

In order to apply compactness results, we need some estimates on Ph,At^h,At an d 
\F h,At■ Let us begin with the estimate on 'J’h.At- 


Proposition 2.1 Under the assumptions of Proposition [TTT1 there exists a constant 
C depending only on the data and independent of At and h, such that: 

ll^ +1 H?,r < C, V &: > 0, (29) 

and ||o;i,T ^ CT. (30) 

The proof is left to the reader. It follows the line of the proof of Proposition 2 
in HU- It mainly uses (|2U1) and (1271) . 

Remark 3 Thanks to estimates (1261) and (l29j) . we have ||||lo=>(o,i) ^ C for all 
k> 0. 

Let us now state the same result on Ph.At and Nh,At- 


Proposition 2.2 Under the assumptions of Proposition [Lll there exists a constant 
C depending only on the data and independent of At and h, such that: 

||«h,At||o ; i,r ^ LI, for u = P,N (31) 

Proof The proof is based on a classical method, already applied in |18j . How¬ 
ever, it requires to pay a special attention on the boundary conditions which induce 
a new difficulty. 

Multiplying (114bl) with A tu k+1 and summing over i and k, we obtain A + B = 0 
with 

A =J2J2 e uhiU^ +1 and B = Atu i +1 

k =0 i =1 k =0 i= 1 


(p^ 1 , -p k+ \ 
V U ’ l +2 U ’ l ~2 


It is easy to see that: 


K -1 I 


k =0 i =1 






i =1 




(32) 
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Moreover, applying a discrete integration by parts to B and using the following 
decomposition of the numerical fluxes given in f[3j: 


■pk+l 


-z u dV k+1 di 


,k +1 
l i +1 




+ 


z cW k+ } 
Zl, 'z+i 


coth 


-zuh lH d* k +r 


, k +} _ „.k+l\ 


H +1 


(33) 


we can rewrite B as B = B\ + B 2 + with 

K ~ 1 1 A tz (-Zuh l+ idd> k +\' 

= ~^ k+ ' coth 1 + 2 ,+T 


i+i 


k =0 i =0 
A-1 I 


-s-oh 1 )' 


A tz,, 


■% (“K) 2 -k +, v 


a > = EE-f d *^ *’ 

k =0 i=0 V 

s 3 = x; At 


/c=0 


As xcoth(x) > 1 for all x E R, we have, as in [4], 

A'— 1 I A , 9 

/ fc+i _ „*+*) . 


B ' > E E 7TT 


(34) 


fc=0 2=0 2 

Applying a discrete integration by parts and using the scheme (HH), we get: 
K -1 1 


a = - V ( ' -* Of 1 )’ 

At "” („s+‘) ; 


fc=0 »=1 
A—1 

+ E 


fc =0 


Since zpzjv (^^ +1 — u max ^j (u k+1 ^J ^ 0 for u = N, P, we have: 

K A A<*A (^ +1 - pm ") + w (ivf +1 - JV 


EE 

fc=0 t= 1 


A 2 


E f7/(=+ 1 ^ 


A—1 7 


> 21E ^4^ (“h 1 - “™“l f«? +I '' 2 

fc=0 i=l 




which yields 


A —1 1 


/j 2 ^ E E Mz " ,hl 


k =0 i=l 


2A 2 



(35) 
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with 


B a = 




Using the boundary conditions (1161) . we may now rewrite B 3 + B 4 as 


K -1 


B 3 + B 4 = -Y, At «fu) k+1 + (fu) k+1 ), 


k =0 


with 


■°i fc + i = L^y 


Uv )' v ' " = l« 


-A 


,fc+A , £« ^o +1 ~ A 4'\ 


° ; 


+ 


/r 


ao 


+ (V>o +1 ) 


9 r t/ _ / fc+i _ a / p 2c i 

(/v +i = (»?h) -a (v - tft!) - 1 ^ + „{+i 7 i (v _ 

It remains to find an upper bound of (/„) fc+1 + (/u) fc+1 - To this end, we proceed as 
in [1] and HD- We introduce: 

£°(z) = 720*0 - u max P° u (x) + u max ^ (x - AC), Vrel, 

ao 

£(*) = 7i(*) - « max /^) - « max - (® - A^f'=), feel, 

which are nonpositive functions under the hypotheses (I18D - Q19I) (see [T]). As 


■0\fc+l _ 


u 


fc+iV 


2 u ma 


■0 (j,k+l\ _ 


C(i- 


rj --’“/s: a + 1 ) -7(tS +1 )l +7(v'o‘ +1 )“5 +1 


clearly have: ( fu) k+1 < 7° (V’o +1 ) n o +1 - Rewriting ( fl) k+1 with the help of , 
similarly prove: (/„) fc+1 < (v — ^ 7 + 1 1 ^ 7 + 1 - Then, using the estimates (1201) 


(/, 

we 

we 

and (f29l) and the continuity of the functions 7 ° and 7 ^, we have 

b 3 + b 4 > — u. 

From (1221), (1271) . (1221) and (1221) . we deduce that 

. 2 


‘St 1 - “? +1 ) 


K -1 I 


K—l I 

EE A1 fc . , 

fc=0 i=0 2 fc=0 i=l 


+ 7 E& 


u . fc+1 - 




(36) 


(37) 


with C depending only on p max ; I\r ma:E , ao> aq, U, A\Hq ZC , A^ 2C , A, (/3„, 7 u ) ?;=0 x 
and T. This ends the proof of Proposition 12.21 □ 


Remark 4 Note that a direct consequence of (1271) is 

K-i 1 2 

£u^ 2Y1 hi (^ +1 “ u y ^ C ' ( 38 ) 

k =0 i=l 

with C a constant independent of h and At. This inequality will be used in Sec¬ 
tion 13.21 
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3 Compactness results and passage to the limit 

In this section, we prove the Theorem 11.11 Firstly, we establish the convergences 
of ( P m )m and ( N m ) m using a Gallouet-Latche compactness Theorem (Theorem 3.4 
in [23]), which is a discrete counterpart of the Aubin-Simon Lemma. Secondly, we 
prove the convergence of (4' m ) m using a classical Kolmogorov Theorem. Then, we 
show the convergence of the traces on the boundaries following the ideas of [6|. 
Finally, passing to the limit in the scheme, we prove that ( P m ,N m , 'F m ) m tends to 
a solution of (V) in the sense of Definition 11.11 

3.1 Compactness of (P m ) m and (N m ) m 

The proof of compactness being analogous for ( P m )m and ( N m ) m , in all this section 
we use the notation (u m ) m where u can be replaced by P or N. To prove the 
compactness of the sequence we will use Theorem 3.4 in |23j. Therefore, for 

any function Wh,At. € Ht, At, we define a discrete time derivative dt,T w h,At and a 
discrete space derivative d x ^Wh.At- There are piecewise constant function in space 
and time, defined by: 



Let us first establish an estimate on the discrete time derivative needed for the 
convergence proof. 

Proposition 3.1 Under the assumptions of Theorem 11.11 there exists a constant 
C depending only on the data such that: 



(39) 


Proof Let consider Vh m € 'HTm such that ^ 1. By definition, we have: 



Then, using the scheme (|14bjl . a discrete integration by parts and the reformulation 
(1331) of the fluxes, we get 
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with 


Ai = — ^2 \ v i +1 - Vi\ 

£u . „ 
i=0 

I 

\Vi+l - Vi I 




/c+1 . „,&+l 


*4e 


i =0 


-Slcoth 1 1 +I 


<‘ + + l‘ - «‘ + ‘) 


^3 = — 


^ + i** 

a, 2 


1 

+ — 

^7/. 




Using (l20li and (l29j) . we obtain: 

(n i+ i - Vi) z 


1 / / f\Drfc + 1 — vI/H 1 ^ 2 ^ 2 




E 

V 2=0 


h, 


*+3 




i=0 


h. 


i+i 


C 

^ — IKJ|l,T m - 

Eu 


Since x i->- xcoth(x) is a 1-Lipschitz continuous function and is equal to 1 in 0 and 
thanks to Cauchy-Schwarz inequality, we obtain: 


1 V—I Vi 1 1 — Vj 

a 2 <—22 


L i. h -1 2 . 
t =0 2 


z v h T+ uW k +l 

J +2 1+2 


-zuh^tr 

coth | - 4 -- I — 1 


u k+1 — vh 

u i +1 


1 |Ui+l - U;| 


_i_ _ ^ Iii 


i =0 ®+ 


< 1 — ii^i +1 i 


u k+1 — a* 

u i+l u i 


l.Tm + — ll M m l|l,T m IKJ|l,T m - 

J 

Let us now consider the term A 3 containing the boundary conditions. Using 
and the continuity of the functions (/?*, r yl l )i=o, 1 , we have 




£ u A 3 sC ( 7 ° ( 

+ (ii ( 


+ 


V-V 


,k +1 


/c +1 
1 +1 




M 


,fc+l 
l I +1 


Pi 


V- T 


k+1 

1 +1 


)) |v/+i| 


<C(l + |a£ +1 | + 


a 


fc+i 

7+1 


(M + |u/_|_i|) 


^ C(l + ||«m +1 lli,r m )ll^ m ||i,r m - 

Therefore, we obtain: 






CIK 


m 111? 7m 


(i + lhm +1 Hl,T m ) 


which yields 
K -1 


E At -»KV„ 


u. 


’ m ll— 1,2, 7rn 




2C* 2 T 2C 2 


+ 


k=0 


E At 

k =0 


.... V+l 1|2 
m ll “m 111, 7^ 


Finally, the estimate (1391) is a consequence of the L 2 (0, T, iL 1 (0,1)) estimate (f3TT) . □ 
We are now able to prove the following convergence result. 
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Proposition 3.2 Under the assumptions of Theorem 11.11 there exists 
u E L 2 (0, T; H 1 (0,1)) such that, up to a subsequence, 

Um —t u strongly in L 2 (0, T; L 2 (0,1)), when m —» +oo 
9x,Tm u m —*■ weakly in L 2 (0,T;L 2 (0,1)) when m —>• +oo. 

Proof The sequence ('H.T rn ,Atm)m is a sequence of finite-dimensional subspaces 
of T 2 (0,1). Each subspace At m can be endowed either with the norm || • || 1 , 7 ^ 
or with the norm || • ||- 1 , 2 , 7 ^- These two norms verify Lemma 1 2. II and Lemma 12.21 
which correspond to the hypotheses of Theorem 3.4 in [23] . 

Then, thanks to Proposition 12.21 and 13.11 we can apply Theorem 3.4 in [23] : 
up to a subsequence, u m —>• u strongly in L 2 (0,T, L 2 (0,1)). Moreover, Proposition 
12.21 implies that u E L 2 (0, T, H 1 { 0,1)). Finally the weak convergence of d x< j- m u m in 
L 2 (0, T ; L 2 (0,1)) is also a consequence of Proposition 12.21 □ 


3.2 Compactness of (\P m ) m 

The convergence of the sequence (T m ) m will be obtained as a consequence of the 
Kolmogorov compactness Theorem, as it is done for instance in [18 ]. 

Let us first remark that the following estimate on the space translates of '&h.At 
is a consequence of the estimate (f30l) . 


Lemma 3.1 Under the assumptions of Theorem ll.il there exists a constant C such 
that for all rj < h: 

||’Jr,Ai(' + I?) •) — ^ft,At(‘) ') IIL 2 (0,T,L 2 (0,1—r/)) — ^ 

Then, we can also establish an estimate on the time translates of ^h,At- 


Lemma 3.2 Under the assumptions of Theorem ll.il there exists a constant C such 
that for all t < At: 

ll^>At(-) • + r ) — ')IIL 2 (0,T—r,L 2 (0,l)) — ^ T 

Proof Using (I14al) . (|16al) and (|16bD . we have, for i E [l;/], 


-A 2 ( dT fc+2 - dT ^ 1 - dtf K+ ? + dT^' ) = 


-ffc+i 


k+2 


rfc+1 


*+ 


hi (3 (P k+2 - P k+1 ) - 


T £ +2 - ^ +1 = a 0 ( dTi +2 - d ^ +1 ) , 


^k+2 _ $ *ti = _ ai 1 d ^*+* _ 


, k -(-1 


1+1 


1 +1 

Multiplying by 

K -1 1 


-,k+2 


t +-2 


yfc+1 


k-\-2 _ 


k =0 i =1 
K—l I 


I and summing over i and k, we obtain A = B with 


jk -\-2 _ 


-A 2 I d*+? - dtf *+? - dtf + dT^ 1 I fa/ 


b = EE'>' 3 ( p , k+2 - p i 


jk+1 


- (iv fc+2 - N k+1 




k-\~2 _ 


* J 


k =0 »= 1 
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Using the boundary conditions, in a same way as in the proof of Proposition ^. II we 
get: 


.4»AV£ii*j+ 2 -<i'!;+ 1 ii5,r, 

k=0 

with a = min(l, l/ao, 1/ai). Then, using Young’s inequality, (1271) and Remark 01 
we also have: 

B «^EK + 2 -'i'!; +1 iii,T+c. 

k=0 

Then, we obtain 

EK + 2-*E +1 || Ir^C, 

k=0 

with C a constant independent of r and At. With (I27|) . this yields the expected 
result, since 


K -1 


\*hM-i -+r)~ Oil V(0,T-t,V{0,1)) = T E II ~ * 


fc+l 112 

h Ho- 


fc =0 


□ 

Finally, we deduce from Lemmas 13.11 and 13 . 21 the following convergence result for 
the sequence ('I f m )m- 


Proposition 3.3 Under the assumptions of Theorem EH there exists 
T G L 2 (0, T; H 1 (0,1)) such that, up to a subsequence, 

—>• T strongly in P 2 ( 0 , T; P 2 ( 0 , 1 )), when m -> +oo, 

9x,Tm^m dx'Zf weakly in P 2 (0, T; P 2 (0,1)), when m -> +oo. 


3.3 Convergence of the traces 

Up to now, we have proved the existence of P, N and T belonging to L 2 (0, T ; H 1 (0,1)), 
such that, up to a subsequence, (P m ) m , (iV m ) m and (T m ) m converge respectively to 
P, N and T. It remains to prove that (P, N, T) is a solution to the corrosion model 
in the sense of Definition 11.11 Therefore, we will pass to the limit in the scheme, see 
Section EH But, at this stage, we will have to pass to the limit in some boundary 
terms. Therefore, we need the convergence of different traces. 

For w = N,P, or T, we have established that w € P 2 (0, T; H l (0,1)). Due to the 
compact embedding of P 1 (0,1) into C([ 0,1]) in ID, it is clear that the trace of w in 
0 and 1 is equal respectively to rc(0, •) and w(\, •), which belong to U 2 (0, T). 

For a given function R’ft.At € Pt.ai, the usual trace, denoted by 7 Wh,At is a 
piecewise constant function of time defined by 

TWh,At(0,t) = w\ +1 and jw ht At(l, t) = w k T +1 , Vt € [t k , t k+1 ). 
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But, it will be easier to deal with an approximate trace ^Wh,At defined by 

7W/i,At(0,i) = ^o +1 and r»7i,At(M) = Wj+l, Vt € [t k ,t k+1 ). 

Proposition 3.4 For w = N,P, or T, the sequence -)) m (resp. ■)) m ) 

converges towards tn( 0 , •) (resp. u>(l,-)) strongly in L 1 ( 0 ,T), up to a subsequence. 

Proof To prove this Proposition we follow the ideas presented in the proof of 
lemma 4.8 of [ 6 j. More precisely , for all g > 0, we write: 


1 

1 


e T 

\jw m (0,t) - w(0,t)\dt = ^ [ [ \^w m (0,t) - w(0,t)\dtdy ^Ti + T 2 + T 3 , 


o o 


(40) 


with 


e T 


Ti = 


T 2 = 


0 0 
f? T 


0 0 
Q T 


T 3 = - 
6 


| 7 t/> m ( 0 ,i) - w m (y,t )| dtdy, 

\w m (y,t) - w{y,t )| dtdy, 

\w(y,t) - w(0, t)| dtdy. 


o o 


Let us first note that T 3 — >• 0 when g — >• 0, as u;(0, t) is the trace on x = 0 of w(-,t). 
Let us now show that T\ anc. 

e 


Let us define K e = 


hr, 


T 2 converge also to 0. 
+ 1. We can show: 


1 K -1 K e 

T\ ^ - ^2 At m ^2 hp 


w k+1 — w k+l 

w 0 w p 


K-l 


K e P-1 


^ a t. m y^h p y~] 


vUi 1 - +1 


k =0 p=l k =0 

Using Cauchy-Schwarz inequality and (1301) or m 


p= 1 i= 0 


T\ y cVf 


{hm + g) 2 


hm —^0 


CVfg 3. 


Moreover, 


e t 


T 2 ^—\J J (w m (y,t)-w(y,t)) 2 dtdy 

Vo 0 

which gives with (140)) . (1411) . (1421) and p tends to 0 , 


hm —^0 


0 , 


(41) 


(42) 


1 

lim / | 7 iu m ( 0 , t) — iu( 0 , t)| dt = 0 . 
hm^f o 7 


This concludes the proof. 


□ 
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Corollary 3.1 Up to a subsequence, the sequence ( r jw m (0, -)) m (resp. (jw m (l, ■)) m ) 
converges towards w( 0, •) (resp. w( 1, •)) strongly in L 2 ( 0, T), for w = P, N or T. 

Proof The sequence (7tc m (0,-)) m (resp. (p/w m (l, -)) m ) converges almost every¬ 
where on (0, T). Moreover, jw m ( 0, •) (resp. •)) is uniformly bounded thanks 

to the Propositions ll.il and 12. 11 Then, we obtain the strong convergence in L 2 ( 0, T). 

□ 


3.4 Passage to the limit 

We end the proof of Theorem 11.11 Indeed, it remains to prove that the limit P, N. 
T, defined in Propositions 13.21 and 13.31 is a solution of (V) in the sense of Definition 
11.11 To this end, we follow the method used in [4] and m- 

First of all we prove that ( P, N, T) satisfy (fT2j) . For u = P,N and ip € D([0,1] x 
[0, T[), we define 


A 10 (m) 
A 20 (m) 
A 30 (m) 
A 40 (m) 

£A(m) ~- 


£q,Ur 


0 J 0 


>0 Jo 
r T r 1 


(x,t) (dt<p(x,t)) dxdt + J e u u m (x, 0)ip(x, 0) dx 
(d X} T m u m (x,t)) (d x <p(x,t)) dxdt, 

n z u u m {x,t) (d x j- m ^ m (x,t)) (d x <p(x,t)) dxdt, 

: [ -^m(l,t))^U m (l,t) -nfliV -J^ m (l,t))) ip(l, t) dt, 

Jo 

+ / (/^(7^m(0,t))yu m (0,t) - 7°(7T m (0,t))) <p(0,t) dt, 

Jo 

- A 20 (m) - A 30 (m) - A m (m). 


Due to Propositions 13.21 13.31 and Corollary 13.11 we have: 

SA(m) —> / / e u u(d t ip ) dxdt + 

m^+00 J 0 J 0 

- [ {Pl( v ~ -7«(^- ^(M))) <MM) dt 

Jo 

~ [ (/5°(^(0,t))R(0,t)-7°(T(0,t))) ^i(0,t) dt. 

Jo 

Thus, it remains to prove that £ 4(7 n) —>-0, when m -A- + 00 . To this end, we multiply 
the scheme (| 14b |) by A tip*, with p\ = p(xi,t k ), and sum over i and k. Using the 
fluxes decomposition (l33ll and the boundary conditions (1161) . we obtain: 


T r l 


e u u(x , 0)<£>(x, 0) dx + 


Ju (d x <p) dxdt 


0 Jo 


Ai(m) + A 2 (m) + A 3 (m) + A±{m) = 0, 
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with 


K -1 7 


M{m) = Y Y £uhi \ u i +1 - u i) Vi 


k =0 i =1 

K ~ 1 1 rt k+1 r x i+1/2 


A J rt 1 r x i+ 1/2 J r 

Y X] / / dxdt - £u^ / 

7—n 1 Jt k •'Xi-U‘2 Jx. 


'*+ 1/2 


— &U 

k =0 i=l Jtk Jx i- 1/ 2 

A'-1 7 

-42 (m) = - X] X] -2— 1 COth 

fc=0 i=0 


r4V(a7,0) da;, 


i=l " •1'*— 1/2 


-^ +i d^+r 

*+2 1+2 


n fc + 1 _ t,^ 1 

U i+1 U i 


Vi+1-Vi) > 


K -1 7 


“?+l +«? +1 


fc=0 1=0 


-4.3 M = X] X] A *m*udtf*+| 2 

^fc+l j_ ^fc+l 

EE 


V?+l - if?) , 


= ^7, 


A — 1 7 fc+1 , fc+1 r t k+1 

U i+l + “i (Iv pfc+l / 

2 *+| Lk 

k =o 1=0 2 7 4 


Xj + I,r) - (p(xi,t k )) dt 


K-l 


Mm) = Y Af ™ [(Puiv - ^MXl - - <0) Vl +1 

k=0 

+ (/92(®g +1 )«g +1 -72K +1 )) ^o 

Following the standard method used in [4], we can prove that 


\Ai-Aio\ —> 0, for i €[1,3]]. 

m—>•+oo 


The only difference with the standard method concerns the terms related to the 
boundary conditions. 

Then, let us compare A 4 with A 40 . We can rewrite 4 4 q under the form: 


140 


K— 1 £fc+l 

M = Y (f j u( V - ^/+i)«/S - 7 i(V - TXl)) ¥>(1 ,t) dt 

k =0 Jtk 


Y (/ 3 °(< +1 )4 +1 -7°(< +1 )) I <p(P,t)dt 

k =0 Jfk 


and, using the continuity of the functions (/?„, 7^)i=o, 1 , the regularity of p, (1201) and 
(|26j) . we obtain: 


|-4 4 (m) - -4 40 (m)| < C'|MIc2((o,t)x[o,i]) m= ^ oo 0. 

Finally, ej^ijn) —> 0 and P,N, \F satisfy (1121) . In the same way, we can prove 

m—>■+oo 

that P, N and *F satisfy (1131) . This concludes the proof of Theorem 11.11 


4 Numerical experiments 

In this Section, we present some numerical results for a test case close to the real 
case given in [T]. Indeed, we have modified some data of the real case in order to 
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have a simplified model which is close to the DPCM model. For the DPCM and 
simplified models, it was performed a scaling relative to the characteristic time of 
the cations. This gives a very small value for the coefficient e. Therefore, even if 
we have proved the convergence of the scheme only for e > 0, we want to study 
numerically the behavior of the scheme for different values of e. More precisely, we 
are interested in its behaviour in the limit e —)• 0. 

4.1 Presentation of the test case 

The test case, introduced here, is inspired by the real case given in [lj. Starting from 
the data given in [lj, we have built a test case which is adapted for the simplified 
model. We only consider a potentiostatic case, which means that V is an applied 
potential. For the different experiments, we will consider the numerical values for 
the dimensionless simplified model given in Table [I] 

Table 1: Dimensionless parameters. 


A 2 

OL 0 

ai 

pm 

N m 

iv p 

;.0 

k N 

kp 

k 1 

k N 

1.1 10~ 3 

0.177 

0.089 

2 

1 

10 8 

10-18 

10 4i 

26.8 

rrip 

rri % 

nip 

m N 

ap 

a N 

Op 

°N 

CL p 

0 

1.45 10" 24 

10 s 

26.8 

0.5 

0.5 

0.5 

0.5 

0.5 

a N 

b p 

b N 


ATf c 

V 

Phi 



0.5 

0.5 

0.5 

-0.866 

0 

0.5 

-5 




Let us remark that for such a choice of coefficients all the hypotheses (d6j), ((TJ) , 
(f8|) . ceo and a\ positive, (fl8l) and (fl9l) [) are satisfied except for the right inequality 
in (USD. Nevertheless, we observe that numerically (12011 still holds. 

As mentioned in the introduction, a scaling relative to the charateristic time gives 
a very small coefficient e, which is the quotient of the mobilities of the densities. In 
practice, e is close to 10 -14 . Then, for the experiments, we will consider different 
small values of e, e S { 0 , 10 - 6 , 10 ~ 4 , 10 — 2 } in order to verify that the scheme is 
asymptotic preserving in the limit e tends to zero. 

4.2 Numerical results 

First, we want to illustrate numerically the convergence of the scheme, for different 
values of e and its asymptotic preserving behavior in the limit e tends to zero. Since 
the exact solution of this problem is not available, we compute a reference solution 
on a uniform mesh made of 4000 cells with a time step At = 10~ 6 for each values 
of e. For the simplified model, there is a stationary state in long-time. Thus, it is 
important to take a small time T in order to compute a solution different from the 
stationary state. Then, for Figures [H [2] and [3] the L 2 -norms in space are computed 
at the final time T = 10“ 3 . 

Figured] shows the L 2 -convergence rate in space of the scheme for different values 
of e. We observe that the convergence rate is independent of e and the scheme has 


International Journal on Finite Volumes 


21 














































Finite volume scheme for a corrosion model 


an order 2 in space, even for e = 0. It is due to the choice of Scharfetter-Gummel 
fluxes for the numerical approximation of the convection-diffusion fluxes. Figure [2] 
shows the L 2 -convergence rate in time for different values of e. We observe that this 
convergence rate is also independent of e and the scheme has an order 1 in time as 
expected. 

In Figure [3] we present, for different values of the time step, the L 2 error in space 
at the final time T = 10~ 3 with respect to e. The asymptotic preserving behavior 
of the scheme, in the limit e tends to zero, clearly appears. 


(a) L 2 error on the cations density (b) L 2 error on the electrons density 



Number of cells 



L 


(c) L 2 error on the electric potential 



Figure 1: L 2 errors with respect to the space step on the densities and the electric 
potential for different values of e and at final time T = 10” 3 . 

In Figure 01 we present the profiles of the two densities and the electric potential, 
in the case e = 0. In order to reach the stationary state, we use as final time T = 1. 
The solution is computed on a uniform mesh with 2000 cells in space and with a 
time step At = 10 -3 . We obtain the expected behavior for the densities and the 
electric potential. 

5 Conclusion 

In this paper, we prove the convergence of a numerical scheme consisting in an 
implicit Euler scheme in time and a Scharfetter-Gummel finite volume scheme in 
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(a) L 2 error on the electrons density 


(b) L 2 error on the cations density 




Figure 2: L 2 errors with respect to the time step on the densities and the electric 
potential for different values of e and at final time T = 10~ 3 . 


space for a corrosion model. The convergence proof is only valuable for e > 0. But, 
as mentioned in the introduction, the dimensionless parameter e is the ratio of the 
mobility coefficients for cations and electrons and it is therefore very small. It should 
be set to 0 in the practical applications. As seen in the numerical experiments, the 
scheme is robust with respect to the value of e. It seems to be asymptotic preserving 
in the limit e —>• 0 (Figures 11121 and l3l). In practice, we can set e = 0 and obtain the 
expected behavior for the densities and the potential (Figured]). 
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